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• Abstract 
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' A numerical method for integrating the equations describing a dynamically coupled 

^ . system made of a fluid and cosmic-rays is developed. In smooth flows the effect of 

I CR pressure is accounted for by modification of the characteristic equations and 

the energy exchange between cosmic-rays and the fluid, due to diffusive processes 
in configuration and momentum space, is modeled with a flux conserving method. 
Provided the shock acceleration efficiency as a function of the upstream conditions 
and shock Mach number, we show that the Riemann solver can be modified to take 
into account the cosmic-ray mediation without having to resolve the cosmic-ray 
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Q^. induced substructure. Shocks are advanced with Glimm's method which preserves 

^ I their discontinuous character without any smearing, thus allowing to maintain self- 

' consistency in the shock solutions. In smooth flows either Glimm's or a higher order 

c/3 . Godunov's method can be applied, with the latter producing better results when 

' approximations are introduced in the Riemann solver. 
> . 
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1 Introduction 



We wish to formulate a numerical method to solve a system of equations 
characterizing a fluid that is dynamically coupled to suprathermal particles 
through the exchange of momentum and energy. Such conditions occur com- 
monly in astrophysical plasmas. The fluid system is an ordinary nonrelativistic 
gas described by the following modified equations of hydrodynamics : 
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where {p,u, Pg,eg) indicate the gas density, velocity, pressure and specific en- 
ergy respectively; i, d index the spatial components and summation over re- 
peated indexes is assumed; 5id is the Kronecker's delta. The gas total spe- 
cific energy is, Cg = + eth, and a 7-law equation of state is assumed so 
that the gas specific internal energy is related to the gas pressure through 
eth = Pg/pilg ~ !)• The inhomogeneous terms proportional to dPc/dx on the 
right hand side of Eq. (EJ-Q account for the effects of the suprathermal pres- 
sure. S is a source term describing the transfer of energy between the fluid and 
the suprathermal component. This may be due to, e.g., particle acceleration 
processes at the expenses of the fluid energy or, conversely, energy losses from 
the suprathermal particles that end up heating the fluid. 

As for the suprathermal component we consider cosmic-ray (heretofore CR) 
particles described by a distribution function, /(x, p, t), which depends upon a 
spatial, a momentum and a temporal coordinate. In what follows p is in units 
of ^mcc\ with rric the CR particle mass, and the normalization of / is such 
that the number density of particles with momentum between p and p + dp is 
dric = Anp"^ f dp. In addition, / is assumed to be isotropic in momentum space 
and evolves according to the following diffusion-convection equation 
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The second and third term on the left hand side of the above equation rep- 
resent, respectively, spatial advection and diffusion with a coefficient, k(x,]9). 
The first term on the right hand side accounts for adiabatic effects and, 
bi{p) = —{dp/dt)ioss, describes the particle momentum change due to energy 
losses associated with mechanical and radiative processes. In addition, Dp{p) 
is the diffusion coefficient in momentum space. The CR pressure in Eq. ([2]) 
and Q is then defined through the distribution function, /, as 

Pe(x) = — m^c" / / /(x,p) ip' + dp, (5) 

O Pmin 

where PmimPmax are the minimum and maximum CR momenta, respectively. 
More specifically, the momentum Pmin marks (somewhat loosely) the transition 
between the thermal and nonthermal components and Pmax is the maximum 
momentum the particles can achieve and still be confined inside the system. 

The CR energy and adiabatic index are given by 
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E, = ATim,c^ / pV(x,p)[(p' + 1)2 (6) 



7c = l + ^. (7) 
The evolution of the CR energy is obtained from Eq. (jlj) and reads 



at 



dx 



—AnmrC 



Pinax 





bif + Dp 



dp J (p2 + i) 
d\nf 




Pmax 



The first hne of Eq. (IHl) refers to the effects of advection, diffusion (with an 
energy averaged diffusion coefficient (k)) and adiabatic compression, and the 
second hne to energy losses/gains introduced above. The surface terms on 
the third hne describe changes in the CR energy due to the flux of particles 
across the low and high boundaries in momentum space. The first of these 
two surface terms is typically negligible whereas the second is important in 
case of efficient shock acceleration and can cause significant energy losses in 
the system. Finally the source term in Eq. ([3]) is related to the change in the 
CR distribution function due to flux in momentum space as 



rpmax / df\ 

E(x) = —inrricC / bme f + D„— ^ dp 

^ ' Jp^.^ \ 'dp] (p2 + i)| 
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where bme{p) now includes mechanical losses only (i.e. radiative losses are 
excluded) . 

The spatial diffusion coefficient introduces a physical scale characterizing the 
particles mean free path due to diffusion, i.e. Xmfp ~ 'i{p)/v{p), where v{p) ~ 
c, the velocity of a particle of momentum p, is of order the speed of light. 
In the following we shall distinguish two different regimes of application of 
the equations namely smooth flows and shock waves. The reason for 

doing so is that for astrophysical systems, Xmfp ^ Xgystem, and the entire 
dynamic range of scales cannot be resolved with currently available computers. 
However, while on the scales that can be resolved by simulations diffusion in 
smooth flows can be safely assumed to become either slow or negligible, this 
is not the case for shocks. 
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In smooth flows (and on large enough scales, A ^ ^mfp) the presence of CRs 
enhances the propagation speed of sound waves but simultaneously causes a 
damping of their amplitude due to CR diffusion [32I31J- In addition energy is 
exchanged non adiabatically between the thermal and nonthermal components 
according to Eq. iQ. These effects arise from diffusive processes and as long 
as the relevant transport coefficients, k and Dp, are defined correctly, they can 
be properly modeled numerically with schemes available in the literature. 

Around shocks, however, the situation is more complicated because the diffu- 
sion process gives rise to an efficient mechanism for transferring energy from 
the flow to the particles. This topic is discussed in detail in several review 
articles p[5] . Here we emphasize two basic and related points relevant for the 
present discussion. Firstly, the dissipation of energy into CRs changes the value 
of the total pressure generated by the shock dissipation mechanism due to the 
different thermodynamic properties of gas and CRs. In addition, as illustrated 
above (cf. Eq. [8]) the escape of high energy particles upstream of the shock 
allows for energy to be removed from the system. This can reduce the pressure 
support in the downstream region, allowing for compression ratios higher than 
the hydrodynamic limit. Finally, the CR pressure gradient produced by the 
CR particles diffusing upstream decelerates the flow approaching the shock 
front. As a result the velocity structure is not a sharp transition anymore but 
is broadened up to scales of order the diffusive scale length of the most ener- 
getic CR particles. This is X^^ipmax) = K{pmax)/u shock, where Ushock is the shock 
speed. This effect creates the so called shock precursor where the upstream 
gas is adiabatically compressed before being shocked. Thus, even though in 
a numerical calculation the CRs do not diffuse out of the resolution element 
during a timestep, shock acceleration can modify signiflcantly the shock jump 
conditions with respect to the simple fluid case (see Section [3]). 

There are, therefore, at least two different limits of interest for solving the 
equations in the case of shocks. One which focuses on the study of 

the diffusive shock acceleration process itself. In this case one requires: (i) a 
kinetic approach in which the evolution of the distribution function in momen- 
tum space given by Eq. (jl]) is followed accurately [2|23] : (ii) enough spatial 
resolution to properly resolve the full range of relevant scales that enter the 
problem, from the thickness of the shock, I shock, to the diffusive scale length of 
the highest energy CR particles, XKipmax)- A number of codes, with different 
levels of sophistication, employing various numerical methods and devoted to 
this type of approach have been developed (e.g. [T3l[T|19|llllHll4|22|16j ). Some- 
time and to various extents they also include the processes that regulate the 
diffusive properties of the medium (e.g. wave ampliflcation and damping), a 
key factor for the efficiency of the acceleration mechanism. By necessity, they 
focus on a very narrow region around the shock, of order Xnipmax)- Com- 
plemented by analytic studies p^25f23p^ . among their ultimate goals is to 
investigate, as a function of the upstream gas conditions, U~ , and the shock 
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Mach number, Ai, the downstream CR distribution function, f~^{p), and the 
efficiency of the shock acceleration process, 77. Here and in the following this 
quantity is defined as the fraction of total pressure upstream of the shock that 
is converted into downstream CR pressure, P^: 

where is related to through Eq. 

There is then the opposite limit to the one described above, which focuses 
on the dynamical effects of CRs in smooth flows and shocks, for systems of 
size Xsystem ^ ^niPmax)- We cau refer to it as the astrophysical limit. This 
approach is more application oriented |17ll31f3Uf28ll29f36|[T8j and does not aim 
at studying f~^{p; U' or ri{U^ In this paper we attempt to design a 
numerical method that serves this purpose, without having to resolve scales of 
order i shock- In doing so, we seek to eliminate all but the essential information 
about the CR distribution function in momentum space, so that a fluid-like 
description is approached. Information at the kinetic level must, however, be 
preserved in two parts of the formulation: (a) when computing shock solutions 
and (b) when computing the time-evolution of the CR pressure. The ffist 
requirement is set because, as already pointed out, a correct shock solution 
can only be obtained with a fully kinetic approach [2][23] . This means that the 
only way to meaningfully include the effects of CRs acceleration in a fluid-like 
approach is to assume that f~^{p; U~ and the shock acceleration efficiency, 
ri{U^ ,A4), are provided independently (e.g. from kinetic models) as part of the 
input. When the acceleration efficiency is high and the contribution from the 
highest energy particles dominate the CR energy and pressure, the flux of CR 
particles across Pmax must also be specified [26] . The second requirement stems 
from the fact that dPc/dt is basically the result of energy losses/gains and 
diffusion of the CR particles. Since these processes are strongly momentum 
dependent, dPc/dt is bound to be different for different shapes of the CR 
distribution function and one ought to be able to account for this. Note that 
the two-fluid approximation alone, in which Eq. (jl]) is integrated in momentum 
space to derive an equation for the time evolution of Pc, is not sufficient for 
these two purposes. 

In order to properly evolve Pc, we divide momentum space in a set of Np 
(~ 10) coarse kinetic volumes or momentum bins and integrate Eq. (jlj) within 
the boundaries of each of them. This provides an equation for the evolution 
of the number density of CRs within each bin and is a cost-effective way of 
following the change of shape of the CR distribution function resulting from 
the momentum dependent CR processes mentioned above. This essentially 
works because the CR distribution, /, is typically a smooth power law with a 
slowly varying slope as a function of momentum. This fact provides a natural 
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way for describing f{p) within each bin [27p7] . which compensates for the 
coarseness of the discretization of momentum space. 

In addition, we show that once the shock acceleration efficiency, 77, is specified 
it is possible to account for the modifications induced by the CRs on the 
hydrodynamic shock solution, even though the structure of the shocks remains 
completely unresolved. This is achieved by solving a slightly more complicated 
Riemann problem, after proper modification of the definition of the nonlinear 
waves that appear in it. This shock solution can still be thought of in terms 
of a two-fiuid model description [TUI but with the important difference that, 
among the family of admissible solutions [2] , we select the one demanded by the 
(explicitly) adopted shock acceleration model. This ensures a self consistent 
description of the CR-hydrodynamic system. 

In order for this to work, however, it is essential that the shock discontinuity 
does not spread as a result of numerical dissipation. This is because in general 
the dissipation of CR energy at a shock is a nonlinear function of the full jump 
conditions. Therefore, if the shock is artificially spread over a few zones the 
sum of the CR energy generated at each numerical subshock will in general 
not be the same as that predicted by the model for the full jump conditions. 

A suitable hydrodynamic method for our purpose is the one originally pro- 
posed by Glimm [15]. Introduced as part of a constructive proof of existence 
of solution to nonlinear hyperbolic conservation laws it was turned into an ef- 
fective numerical scheme for hydrodynamics by Chorin piTj . Glimm's method 
maintains unsmoothed all the sharp features that are present in the fiow. In 
particular, shocks remain unsmoothed jumps as they propagate across the 
grid. This allows us to maintain self-consistency in the shock solution. The 
limitation in using Glimm's method is that at the moment its multidimen- 
sional extension does not work properly at shocks [S]. Thus here we focus on 
a one dimensional algorithm and leave its generalization to more than one 
dimension for future work. 

In smooth fiows either Glimm's or Godunov's method can be applied. There- 
fore it is possible to define a hybrid scheme where Glimm's method is applied 
at shocks and Godunov's method in smooth parts of the fiow. In either case, 
(in smooth fiows) the effects due to the CR pressure are included by proper 
modification of the characteristic analysis. 

Note that recently a method has been proposed in [T51I53] to include the dy- 
namical effect of CR pressure on the hydrodynamics. That approach consists 
effectively of a two fluid model in which the generation of CRs at shocks is 
treated as an explicit source term, similar to the scheme in [27]. However, a 
hydro scheme that includes self-consistently the kinetic effects on both the 
shock substructure and the CR pressure evolution, in the sense mentioned in 
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(a) and (b) above, is still lacking. 

This paper is structured as follows. In Section [2] we describe the discretization 
of momentum space and compute the fluxes due to energy losses and diffusion 
in that dimension. In Section [3] we discuss the effects of CRs on the structure 
of the flow and define a modified Riemann problem to include the effects of 
CRs both in smooth flows and at shocks. In Section 14.11 and 14.21 we describe 
the implementations of Glimm's method and of a hybrid Glimm-Godunov's 
method, respectively. Tests follow in Section [5] and conclusions are presented 
in Section [61 



2 Diffusion-Convection Equation 



In order to formulate a finite-kinetic- volume description of the diffusion-convection 
equation (jlj), we divide momentum space into Np logarithmically spaced bins, 
each with boundaries p,- , i . The log-width of the bins, ^Wj = log(n,_, i ), 
is taken as constant (although this is not necessary). We then follow the evo- 
lution of the CR number density associated with each bin, namely 

rip^ = r-'Unp' f{p)dp. (11) 
For a piecewise power-law distribution function we have 



/(p) = /,(p) = /( 



P 



n. 



(12) 
(13) 



where /oj and qj are the normalization and logarithmic slope for / in the jth 
momentum bin. Once the set, [rip.] < j < Np}, is defined and the boundary 
conditions for the slopes g_i and are provided, we can compute the set of 
slopes, {qj] < j < Np}, and normalizations, {foj; < j < Np}, as follows [T7] . 
For each bin, i, we use Eq. f ll3p with j = i,i ±1 and further assume: (a) that 
the spectral curvature, dq/dlnp, remains constant across adjacent bins; (b) 
that f{p), as given in Eq. f|T2l) . is continuous across the bins boundaries, 
and Pj+i . This provides six equations for six variables that can be efficiently 
solved with an iterative method. This procedure is applied for each bin, and 
allows to reconstruct the set {fj{p)} from {up.}. In the following we will use 
both Upj and fj{p) assuming that we can reconstruct the latter from the former 
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through the procedure just describe 



The equation describing the evolution of, Up. , is obtained by integrating Eq. (0 
multiphed by a factor 47rp^, within the interval, — p, , i - This leads to 



^ + V, ■ = -VpFp + J„ (14) 

F^ = u Up^ - (k) Vup. , (15) 

Fp = 47rp^f,{p)p (16) 

where Vp is the undivided (one-dimensional) gradient in momentum space 
and we have introduced 



(k), = M pWfjip) dp I M p'vfjip) dp, (17) 

1 d\n f 

p = - (V-n)p-h{p)-Dp- 



dp 

In writing Eq. f[T^ - f[TB]) we have emphasized that both advection and diffusion 
terms along the momentum coordinate can be cast in conservative form, in 
analogy to the corresponding terms in configuration space. This allows us to 
adopt a Godunov-like scheme for the numerical integration of those terms. 
However, we place them on the right hand side of the Eq. ( fT4l) because they 
will effectively be treated as source terms of Up^. Finally, on the right hand 
side of Eq. ( HM . we have added a source term, Jj, which represents the rate of 
production of CR particles due to the diffusive shock acceleration mechanism. 
We do this only pro forma: because we always use Glimm's method to advance 
shocks, the Riemann solver effectively subsumes the role of the term Jj. In 
other words, Jj will not be treated as an explicit source term but as an implicit 
part of the shock solution computed through the Riemann solver. Thus, it will 
be sufficient for our purposes to only formulate a prescription for the postshock 
values of the set {up.}, which we do in Section [3l3l 



^ When spectral curvature is important, the following alternative approach first pro- 
posed in [27] can be used instead: for each bin, in addition to the number density of 
CR particles, one also follows the energy density, €p.. For each bin the definitions of 
Up. and Cp. provide two equations which can be readily solved for the two unknowns 
{foj, Qj)- While the equation for Up. is derived below (cf Eq. [T^ . the one for ep. is 
obtained analogously by multiplying Eq. ([ID by a factor, 47rp^[(p^ -|- 1)2 — 1], and 



integrating it within the interval, — See also 1161 for the effectiveness of 

2 ■'"'"2 

this method. Here we take the simpler approach, however, in which we follow Up. 
only. This is because we wish to focus on the novelty of the method which is related 
to the fluid aspect of the solutions. In fact, the method for the evolution of f{p) in 
phase space was extensively studied in [27] . 
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2.1 Fluxes in Momentum Space 



Time integration of Eq. (fT4|) due to the fluxes in momentum space is done 
following the method proposed in [27]. We can also retain the diffusive term 
although here we limit the discussion to the case of a small diffusion coefficient, 
Dp, which can be treated explicitly, i.e. 

TD, = ^ » At. (19) 

Up 

Here rn is the characteristic diffusion time in momentum space, A» = , i — 
the momentum bin size, and At the time step. Then, the terms on the 

-' 2 

right-hand-side of Eq. (fT^ produce a change in Hp. such that 



<f ' - < = (d - Ci ) - [(V. ■ F.) + J,], (20) 



^3 •'-^ 



where the superscript n + indicates time centering and the last term on 
the right hand side will be specified in Sections 14.11 and 14. 2[ We can start 
estimating the time-averaged flux in momentum space at time t = nAt as [27] 

^^.^=^r P" flip) dp, (21) 

which is obtained by time integrating Eq. ( fT6l) and by changing integration 
variable from time to momentum [27]. In Eq. ( l2Ti) is the upstream distri- 
bution function at time t defined at the grid point in momentum space 

?■ + 1 if pi.i < 0, 
ju={ (22) 
3 if P^+i > 0, 

and Pu is the upstream momentum, solution of the integral equation 

At=r^. (23) 
^^..^ P 

Note that the denominator of the above integrand function has typically a 
polynomial form so that the integral can be computed in closed form [27J. 

Finally, a time centered estimate of the term, V-Fp ^ , is obtained by taking 
the average of the fluxes, Fp, computed at t and t + At, as usually done 
for nonstiff source terms. Time centering is relevant because the function in 
Eq. (!23|) in general depends on the local properties of the fluid, such as density 
and temperature, which change with time. Note that there are two ways to 
compute the divergence of the velocity in the p term: with a cell centered 
scheme (Vu)i = (wj+i — Mi_i)/2Ax or with a staggered scheme, (Vm)^ = 
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(Wj+i — )/Ax in which case Wj^i is computed as part of the solution to the 
Riemann problem that one has to solve in order to advance the fluid equations 
and estimate the spatial terms appearing in Eq. (1201) . This is described in the 
next section. In the following we always use the staggered scheme except in 
the Godunov's predictor step. 



3 Riemann Problem for Cosmic-ray Hydrodynamics 



In this section we describe the modifications to the Riemann problem in the 
presence of CRs. Without loss of generality we restrict to the one-dimensional 
case. In addition, for the sake of clarity, in the following discussion we shall 
neglect the spatial diffusion term, except for the fact that it is implicitly as- 
sumed to be at work at shocks, where it causes CR particles to be accelerated. 
In smooth flows this term is assumed to be slow and is taken into account 
with an explicit conservative formulation. 



We begin with rewriting our system of equations in conservative form: 

S{U), 



9U ^ ^ 
^ + V, ■ 
ot 



(24) 



where 



U 



pu 
pe 



( 



pu 

pU^ + P 

(pe + P)u 



V 



Up^U 



S 



I 






E, 

-'-'loss 



(25) 



the entry Hp. is repeated for all momentum bins, e.g. for j = 0,Np — 1, and 
we have introduced the total pressure and specific energy as 



P = Pg + Pc 
1 2 



pile, - 1) pile - 1) 



(26) 
(27) 



Losses in the total energy of the system are due to escape of energetic particles 
and, in principle, radiative losses. Thus we write 

Ei^ss = -47rm,c2 p^fp (p^ + 1)2 - I _ + & dp , (28) 

[ L J P-Pmax Jpmin (^^ + 1) 2 J 
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with br referring to radiative losses only. Note that the first term is different 
from zero only if p > because we assume f{p) = for p > Pmax- In addition, 
just hke the term Jj discussed in the previous section, the source term Eioss 
will not be treated explicitly but will be part of the shock solution computed 
through the Riemann solver. 

The conservative character of the Eq. fl24|) for the system in Eq. fl25l) suggests 
that the jump conditions across a shock wave can be written in the usual 
way, provided that the total, thermal plus cosmic-ray, energy and pressure 
are used. However, there is an additional complication related to the way in 
which the total pressure and energy is partitioned between CRs and thermal 
components. This is further addressed in Section 13. 2[ 

For the sake of clarity in the ensuing discussion, we now outline the Riemann 
problem. Suppose that at t = the gas is described by 

U = { - (29) 

[U^ if X < 0. 

The solution at t > is in general characterized by two waves: a backward 
moving wave separating the states {U\ f/*') and a forward moving wave sep- 
arating the states {If^, U*^). Each wave will be either a shock wave or a rar- 
efaction wave. The central states U*\U*'^ are separated by a slip line across 
which the velocity and total pressure (m*, P*) remain constant, but the density 
and individual pressure components, {p*° , P*° , P*°) , o = l,r, in general will 
change. The value of these three quantities in each intermediate state, U*°, 
will be reconstructed from U° and the type of wave connecting the two states. 
In the following two subsections we describe the structure of such waves. 



3.1 Characteristic Analysis 



Transforming into primitive variable space we obtain 



U = 



pu 
pe 



'p ^ 



V = 



u 



Pr 



(30) 



where we have introduced the CR concentrations, yp^ = np.rrip/ p. These quan- 
tities are followed as passive scalars and, for simplicity, will be omitted in the 
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following analysis. Note, however, that the source term ensures that there is 
consistency between their evolution and that of Pc, to which they are related 
through Eq. ([n]), (US]) and (P. 



The system of equations for the primitive variables, V, is obtained with the 
following transformation 



Sv 



A{V) = VuV ■ VuF ■ VvU, Sv = VuV ■ S, 



(31) 
(32) 



where 




(33) 



P 



(34) 



Here Cg and Cc correspond to the speed of sound associated with the gas and 
CR pressure respectively. (In principle, a coefficient {d In Pjd In p)s should be 
used instead of in the definition of Cg, but the difference is negligible for 
the purpose of this paper.) Solving for Det{A — XI) = 0, the eigenvalues of A 
are found to be Aq = u — Cg, Xi = u, X2 = u, X3 = u + Cg- The associated left 
and right eigenvalues are, respectively. 



-p/2cs l/2ci l/2c 

1 


pI2cs 1/2c2 l/2cl I 



(35) 



and 
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R 



1 1 1 ^ 

-Cs/pO Cs/p 



0-1 cl 



(36) 



According to this analysis 'simple waves', including rarefaction waves, are 
described by the following equations 



dp 
dP 
du 



cr 



1 

dP ^C^' 



dPc 
dPr. 



(37) 
(38) 
(39) 



where C# = pc#, ^ = g,c, s, is the Lagrangian sound speed or mass flux 
across the wave. The first two equations fl571) - fl551) are the usual relation for 
hydrodynamics with the sound speed modified to account for the CR pressure. 
The last equation describes the change of CR pressure as a function of the gas 
pressure during an adiabatic process. 



3.2 Cosmic-ray Mediated Numerical Shocks 



In the following we consider the structure of a shock modified by the presence 
of accelerated CR particles. As in the case of hydrodynamics, we assume that 
steady state conditions have been reached. When CR acceleration operates 
this takes of order the timescale to accelerate thermal particles to relativistic 
energies. While substantially longer than for a pure hydrodynamic case, in our 
astrophysical limit this is typically still much shorter than the size of a time 
step (cf. [H] for further details on this). 

If we label quantities far upstream and far downstream of the shock with, — 
and +, respectively, the Rankine-Hugoniot conditions read 
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p+ _ p- 

u+-u- = ± — , (40) 



p+ -P' 




(41) 



-^"-"--2 1/p.-l/p- ' (^2) 

where W is the Lagrangian speed of propagation of nonhnear waves. These 
do not yet specify the amount of CR energy and pressure generated by the 
dissipation mechanism. In fact, Pc and 7c are not and cannot be specified by 
the above equations alone. To do that we use the following facts without proof 
(but see [T0|l3ll2] ). In the presence of CRs the shock discontinuity is replaced 
by a precursor, where the gas is compressed adiabatically by the upstreaming 
CRs and Pg (x (x u~"'] the precursor is immediately followed by a viscous 
subshock where entropy is generated but both the CR pressure and CR energy 
flux remain continuous. This means that the structure of the subshock is purely 
hydrodynamical. If we use the label, 0, to indicate the quantities just prior to 
the subshock, the compression at the precursor is given by 

r, = ^. (43) 



Define the shock Mach number as 



M = ^. (44) 



The gas pressure jump across the total shock transition, being the result of 
the precursor adiabatic compression and subshock compressiorlfj, reads [2] 

2^j2^Af_2^ ,45) 

Pg 79 + 1 '^P 7c, + 1 

An analogous relation can be obtained for the conservation of mass equation, 
namely 

Using these results with Euler equation gives [2] 

^ = ^ + 1 - + 7.(1 - r^')M'. (47) 
9 a 



^ A term describing non adiabatic heating in the precursor due, in particular, to 
Alfven wave dissipation, can in principle also be included in Eq. (|45|) . It is neglected 
here, however, for simplicity. 
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Recalling that, M = W/C, Eq. (jl5]) and (jUj) then lead to the following 
modified definition of the Lagrangian speed of nonlinear waves 

where C~ = {'^gp^ P^)^ , P = Pg + Pc. Finally, the value of the CR adiabatic 
index downstream of the shock, 7+, is determined by the energy equation 
through the following relations: 




7c 



79^+-7a + l' 



Q 



loss 



7g-l -1 , ^_IV1 
7g+l P 79+1 M'^ 



- _ 7c (7g - 1) 
79 (7c" - 1)^ 



Q 



loss 



ElossdXj 



f^ + l-r?+7,(l-r-i)A^2 



(49) 
(50) 

(51) 
(52) 



where Eioss is given in Eq. (I28l) . The term Qioss < describes the energy losses 
occurring at the precursor and shock front. It becomes important when the 
CR acceleration efficiency is very high. In this case, this term must also be 
specified consistently with the acceleration efficiency by the kinetic solution. 



3.3 Riemann Solver Procedure 



Having specified the form of the rarefaction and compression waves modified 
by the CRs, we can now define the procedure for solving the Riemann problem. 
First note that, provided the shock acceleration efficiency, ri{U~ , Ai), as a 
function of the upstream conditions {U~) and the shock Mach number {A4), 
we solve Eq. fH7j) to derive a similar function for the compression at the 
precursor, 

rp = rpiU-,M). (53) 

Given the left and right states in Eq. fl2I?l) . we then want to compute the 
intersection point, {u*,P*), of the two wave curves passing through U\U^ 
in the P — u plane. For this we use the iterative technique proposed in [37], 
with the two shock approximations |8j and additional modifications which we 
describe next. 
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In the absence of shocks we change the Lagrangian speed of nonhnear waves 
as follows in order to account for the CR pressure 



W = C!'^ ( 1 + 



73 + 1 P* - P' 



27g 



(54) 



with Cl'^ = p^'^cY ■ If a shock is present we instead use W given in Eq. fj48l) 
with P"*" replaced by P*. However, unlike the pure hydrodynamic case, W 
now also depends on rp{jvl). Thus, using Jvl = W/Cg, we have the implicit 
equation 



W = W 



(55) 



which also needs to be solved iteratively. In addition the tangent slopes to the 
wave curves in the P — u plane, which are used in the iterative procedure to 
find P*, are modified according to 



dP* 



du* 



2W^ 



2W" 



+ YC 



2' 



Y 



p 



7s + 1- (79- 



(56) 
(57) 



To summarize, the iteration procedure is now given by : 

u = 0, Po* = [CiP^ + CIP' - C\C:{u^ - u^)]/{C\ + CD 
while not converged do 

if P;_i > P''^ then 

wl^^ = wl^^{p;_,,wl'^/cY) 

else 

wt'^ = wt'^iP:_,) 

= 1, r''" = 1, C"'" = ci'"- 
end if 

Z^''' = Z^''-{Wl''^,C^''',Y^''') 

^^i _ ^p*_^ _ piyw^ u*''- = - (P;_i - P''-)/W 
P* = p*_^ - [Z^Z'-/{Z^ + Z'-)] {u*''' - u*'^) 
end while 

u* = {wy + Wlu" + - p^)/{Wl + Wl) 



The criterion for convergence can be set to be, |P* — P*_i\/P* < e, where 
e is a parameter that sets the error tolerance. Note that at the end of the 
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above procedure, in addition to P* and u*, we have also solved for W'"'^ and 



therefore, A^^'^ and Tp'''. 



Once the left and right moving waves have been determined, we proceed as 
follows [8]. In searching for the solution at a given point, ^ = x/t, we set 
a = sign(^ — u*) and define 



, {U\W\c\M\rLy) if (T > 0, 

u° = au°, i = a^, u* = au*. (59) 

We then complete the definition of the intermediate state f/*°. If the latter 
is separated from U° through a rarefaction wave, knowing P*, P° and P° we 
can use Eq. (15^ to estimate P*°, P*". This amounts to solving for P*° the 
nonlinear equation 

/p*o\ 7?/79 

P* = p;° + po (^_|_ j , (60) 

and then setting 

In Eq. fl60|) . 7° is the CR adiabatic index of the U° state which remains 
unchanged during an adiabatic process. The density is then estimated through 
the polytropic law 

Finally, n*° = n°^(p*°/p°). (Note that, in Eq. ([10]) and ([H2D, as in the definition 
of the sound speeds, Eq. the quantity (9 In Pc/<9 In p)^ should be used as 
adiabatic index.) 

In the case of a shock wave, on the other hand, knowing both r° and Ai°, 
we estimate P*° with Eq. P*° with Eq. and p*° with Eq. (jH]). In 

addition, 7*° is defined by Eq. (149|) and, with f~^{p; U",A4°) specified by the 
input kinetic-model, the downstream number density of CRs in each bin is 
given by 



T'^^ Anp^ f+{p; f/°, M°) dp. (63) 

Note that the consistency between ?7(f/°, A4°), f^{p', U°, M.°) and Qioss ensures 
consistency between the values of P^", 7*° and n*° also. 
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We then evaluate c* through Eq. (IMI) and define the wave speeds 



A°, A* 



u° + c°, u* + c* if P* < P°, 



if P* > P°. 



(64) 



If ^ hes ahead or behind the o-wave we can set the solution to 



Pi ^1 Pgi ^ci 



p*°,n*°,p;°,p;°,n;°ife<A*, 

p°,<P;,P°,n°^ ife> A°. 



(65) 



However, if A* < ^ < A°, we have to evaluate the solution inside a rarefaction 
wave. This requires integration of the system (1371) -( |39l) . which cannot be done 
in closed form and can be expensive. An alternative method is to linearly 
interpolate between the states f/° and U*° as 



f/ = (?7*° + (1 -C)f/°, (66) 

C = i-^. (67) 
A°-A* 

This works just fine for Godunov's method. However, we find that for Glimm's 
method, when strong rarefactions are involved, it is important to use the exact 
approach in order to avoid spurious effects. 

Before concluding the description of the Riemann solver we point out the 
presence of a slight inconsistency in the formulation. In fact lim^ ^iW ^ 
Cg, i.e. the speed of weakly nonlinear waves does not tend to the speed of 
sound waves. This is a consequence of the large gap in physical scales between 
the sound waves driven by the total CR and thermal pressure and the CR 
mediated shock waves. The former are in fact long wavelength perturbations 
on which scale the diffusion is slow and unimportant. Such scales are much 
larger than those characterizing the structure of a shock. Thus the conflict in 
trying to reconcile the two solutions is due to the impossibility of following the 
intermediate scales. Since the nonlinear effects due to the process of CR shock 
acceleration are only important for strong shocks, a natural way of solving 
the above conflict is to assume that a shock solution is adopted only if the 
shock propagation speed exceeds that of the sound speed given in Eq. 
This leads to the following condition 

AP > Jl±lJ^p^. (68) 
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4 Implementations 



The method described in the previous section for the solution of the Riemann 
problem can be employed to construct time dependent solutions to the system 
of equations ©-dl]) on a grid. 



4-.1 Glimm's Method 



The first implementation we describe is based on Glimm's method. Follow- 
ing [8], here we simply outline the main procedure that we use and refer the 
reader to the original references p3][6H71l8] for a detailed description. 



Consider a piecewise constant approximate solution at time t" = nAt 



X < {i + ^) Ax, ieV, (69) 



where Ax is the mesh size, At the timestep and V the computational domain. 
We seek to advance the solution by one timestep to t = (n + l)At. To do that 
we solve the Riemann problem at each cell interface, i — |, with left and right 
states given by L^Li and f/f . Denote the solution with 



X 



\)Ax 



t - nAt 



(70) 



where we have made explicit use of its property of self-similarity. If the choice 
of the timestep is sufficiently small, say 



1 , 
At<--—, (71) 

A^ax = max(|<|+c^J,V?Gl), (72) 

the wave solutions of the Riemann problems at each cell interface will not 
interact with each other. Then the set of Riemann solutions, n,^ G T^}, 

each covering a region, [i — l)Ax < x < iAx, defines an exact solution, 
U^'"'{x,t), to the initial value problem in ( 1701) for the time interval, t"' < t < 
fn+i^ The solution at each grid point, i, and time, t""*"^, is obtained by random 
sampling W^'^ as follows: evaluate the solution at the point, x = {i — 1 + 
a^~^^)Ax, within the region covered by 7^j_i „, where, a""*"^, is a randomly 

chosen number, a""*"^ G [0, 1). Note that x E i if a""*"^ > | and x G (i — 1) if 
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a"^^ < i. We can then define 



Glimm 



7^, 



1\ Ax 
At 

l\ Ax 
At 



if a'^+i > |, 
if a'^+i < |. 



(73) 



Following [8J, we use a sampling procedure that is based on van der Corput's 
pseudo-random sequence, so that a"' is the rith element of that sequence. 



4-2 Hybrid Glimm- Godunov's Method 



Shock waves are the only features that need to be propagated without numeri- 
cal smearing. Therefore, we have also implemented a hybrid scheme which uses 
Glimm's method to advance shock fronts and Godunov's method for smooth 
parts of the flow. Having described Glimm's method in the previous section, 
here we briefly outline a scheme based on the higher order Godunov's method. 

In Godunov's method the solution is updated with a conservative scheme 

where the source term has been described in Section O 
The fluxes at the cell faces are given by 

d = P (Ci*) . (75) 

where V. i'^ is obtained by solving the Riemann problem (discussed exten- 

sively in Section [3]) with left and right states (Vi_+, Vi+i^_). These states corre- 
spond to up- wind time averages, which allow to achieve second order accuracy. 
They are reconstructed from the cell center, taking into account the effects 
of spatial gradients and the source term, as follows. At each grid point, i, we 
compute centered and one-sided slopes and use van Leer's limiter to make the 
final choice about the local slope, AVi. Then, the up- wind, time averaged left 
(— ) and right (+) states at cell faces are 

v.,± =vr + l (±/ -^^) T^±im, (76) 

V^,± = V^,± + ^ S^,. (77) 
Here A is given in Eq. / is the identity operator and n indicates the 
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time-step corresponding to time t. In addition 



r±{V)= Y: ilrV)-r, (78) 

±Xj>0 

projects out from the state V the components carried by characteristics that 
propagate away from the cell interface {lj,rj are the left and right eigenstates 
respectively and Xj is the corresponding eigenvalue, described in Section 13. ip . 

In conclusion, our hybrid scheme can be summarized as 

(f/"+^)Giimm if i is shocked, 
{Ui )Godunov otherwise. 



n+l 



We say that the cell, i, is shocked if a shock is going to cross it during the next 
timestep. In order for this to happen at least a shock moving with speed Ug 
with respect to the grid must be present at the interface, i — sign{us)^. The 
criterion for deciding whether or not a wave across an interface qualifies as a 
shock will be based on the strength of the pressure jump across it, \Pi+i — 
Pi\/ min(Pj+i, Pj), and shall take the condition (1^ into account. 



5 Tests 



We now present a few tests illustrating the performance of the methods de- 
scribed in the previous sections. The tests consist of a set of Riemann problems 
with initial conditions 

f {p\ u\ Pi Pi 7^, f\p)) iix< 0.5, 

[(p^u^p;,p;,7e^r(p)) if a; > 0.5, 

for which we compare the numerical and the 'exact' solutions. The 'exact' 
solution is obtained by solving the Riemann problem as outlined in Section [31 
numerically but without any of the approximations involved in the Riemann 
solvers for the numerical methods. In particular, no two shock approximation 
is made, and the exact expression for the speed of rarefaction waves is used. 

In order to allow for an easier comparison with the 'exact' solution we only 
retain the adiabatic terms in the diffusion-convection equation, i.e. we neglect 
the terms Dp in Eq. flTB]) and hf {p) in Eq. flTS]) . Similarly we let momentum 
space range over 15 orders of magnitude with Pmin = 10~^ and Pmax = 10^°. 
This choice, while unrealistic, is made in order to minimize energy losses due 
to fluxes across boundaries in momentum space when studying rarefaction 
waves. 
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For simplicity the initial conditions for the CR distribution function are spec- 
ified as an unbroken power-law, 



with Pmin < P < Pmax- CR particlcs return to the thermal pool for p < Pmin 
and escape the system for p > Pmax- Due to the lack of energy losses the 
evolution of f{p), followed with the scheme presented in Section [2|, becomes 
trivial. However, the accuracy of that scheme has already been extensively 
tested in [27] , so that here we focus solely on the quality of the hydrodynamic 
solutions. 

In solving the Riemann test problems with the numerical methods presented 
in this paper we always employ a grid of 128 mesh points on a domain of size 
unity (so that Ax = 1/128). We use jg = 5/3. In addition we use Np = 16 
momentum bins. Throughout the Section time is expressed in adimensional 
code units. 

As already mentioned, when evaluating the solution inside a rarefaction wave 
with the Riemann solver, we have a choice of either integrating directly the 
Eq. (jnZD-dSnD or use the approximate Eq. (1^ . When using Glimm's method we 
test both approaches and compare the results, whereas when using Godunov's 
method we only employ the approximate approach which turns out sufficiently 
accurate. 

For shock waves, we adopt the following simple prescription defining the shock 
acceleration efficiency 



in which the fraction of total momentum impinging on the shock and dis- 
sipated into CR pressure depends solely on the shock Mach number. In the 
above expression, Aimin and A^s are a threshold and scale parameter, respec- 
tively. While clearly a simplification, the functional form of rj, with a sharp rise 
for A4 > Aimin, followed by a fiattening for > A^^, is partially inspired by 
thermal leakage models and the numerical results described in [20]. We take 
A = 0.8, M.S = 5.77 and M.min = 1-5. In this simplified model we use values 
of Qioss < that, based on the prescribed acceleration efficiency and the en- 
ergy equation, allow for 7c > 4/3. With the above choices, the resulting shock 
solutions always admit a subshock, i.e. completely smooth shock transitions 
do not appear. 

The accelerated CR distribution function is assumed to be an unbroken power- 
law as in Eq. (ISTI) . Using Eq. ([S])-(17]), we can thus write the following relation 



(81) 




(82) 
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Fig. 1. Open symbols: numerical solutions obtained with Glimm's method for two 
shock problems. Solid line: 'exact' solution. Left: solution at t = 0.05 for a left mov- 
ing shock with A4 = 20 (initial conditions in Eq. [85]). Right: solution at t = 0.051 
for a left moving shock with Ai = 5 (initial conditions in Eq. |86j). 

between the slope and the CR adiabatic index 



g = 3(l + a[g][7c-l]), (83) 

Pmin 

in which a{q) > is a function of q through Pc{q) and f{p)- The above relations 
imply that q takes a value in the interval (3, oo) as 7c ranges between 4/3 and 
5/3. 

It should be pointed out that in a realistic shock the slope of the distribution 
function, q, in general depends on p and is determined self-consistently with 
the velocity profile in the precursor and subshock. However, the simplification 
made here about the distribution function, as well as other naive assumptions 
made earlier in this section, are solely for the sake of simplicity or easy compar- 
ison with exact solutions. Nothing prevents the use of more sophisticated and 
realistic kinetic models for //(A^), Pmin, Pmax, q{p) with the numerical method 
presented in this paper. 



a{q) = l 



3PJa) 



{p'fip) [{p' 



+ 1) 



5. 1 Shocks 



We begin with three shock problems. Two shocks with mild upstream ratio 
of CR to thermal pressure, one moving leftward with Ai = 20 and the other 
moving rightward with Ai = 5. And a shock with upstream ratio of CR 
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to thermal pressure equal to one, and moving leftward with Ai = 10. The 
problems are specified by the following parameters: 



Qioss = -155.3796, Qioss/F-{pe) = 0.461045724 

= 1.0, = 20, = 1.0, 

= 0.3, 7^ = 1.34, (85) 

= 12.8305315, M*" = -3.80751023, P; = 103.280692, 

PJ = 512.726578, = 1.33433, 

for the first problem, 

g^os. = 0.0, Qioss/F-{pe) = 0.0 

p' = 4.53983646, = 5.03312097, Pj = 18.1559788, 

= 15.6326773, 7^ = 1.34218, (86) 
= 1.0, = 0.0, P; = 1.0, 

P; = 0.3, Y = 1.34. 

for the second problem, and 

Qioss = -19.93444, Qioss/F-{pe) = 0.222258137 

p' = 1.0, = 10, P^ = 1.0, 



Pj = 1.0, Yc = 1-35, (87) 

p*^ = 7.64274954, u'' = -1.22076909, P; = 42.8536158, 

P; = 104.00589, 7; = 1.33433, 

for the third problem. The initial CR distribution functions are specified by 
Eq. (ISTI) . with slope g'''' determined respectively by 7'''' through Eq. (1831) . In 
addition to the initial left/right states of the Riemann problems we have also 
specified the parameter Qioss defined in Eq. (jS^D and its ratio to the upstream 
energy fiux, F~{pe). So, in the above three examples about 46%, 0% and 
22%, respectively, of the energy flux through the shock front is carried away 
by escaping CR particles with p> p 

max ■ 

In these tests the role of Godunov's method in the hybrid formulation is trivial. 
Therefore we only show the results obtained with Glimm's method. The left 
and right moving shocks described by Eq. (l85l) - (!86l) are presented in the left 
and right panels of Fig.[Tl respectively, while Fig. [prefers to the case described 
by Eq. (!87|) . All plots correspond to a solution time, t = 0.05. For each plot 
the four panels show, from top to bottom, gas density, velocity, gas pressure 
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Fig. 2. Numerical solutions obtained with Glimm's method (open symbols) and 
'exact' solution (solid line) at t = 0.05 for for a left moving shock with = 10 and 
upstream CR to thermal pressure ratio equal to one (initial conditions in Eq. |87j). 



and CR pressure. The numerical solution reproduces the 'exact' solution very 
well, without oscillations or artifacts, despite the fact that the CR pressure is 
comparable or significantly higher than the thermal pressure. Note that the 
front of the right moving shock is displaced with respect to the 'exact' solution 
by one cell. This is a characteristic of Glimm's method: as it advances the shock 
front in discrete steps of size Ax, it may inevitably place the shock on a grid 
position that is offset with respect to the 'true' position. By using van der 
Corput sequence, however, the offset is at most one zone and it eventually 
becomes negligible when compared to the distance traveled by the shock [8]. 



5.2 Rarefactions 



We now turn to the following initial value problem: 

= 0.251188643, = -2.01959396, = 0.1, 
= 0.15703628, 7' = 1.34, 

(88) 

p' = 1.0, u' = 0.0, p; = 1.0, 

pj = 1.0, = 1.34, 
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Fig. 3. Open symbols: rarefaction wave solutions at, t = 0.1504683, obtained with 
Glimm's method using the approximate (left) and exact (right) rarefaction solution 
in the Riemann solver. Solid line: The solid lines indicates the exact solutions. The 
initial conditions are described in Eq. ()88p and correspond to a rarefaction wave in 
the A+ characteristic family. 

representing a rarefaction wave in the A+ characteristic family. The initial 
CR distribution functions are power- laws, Eq. (IHTl) . with the same slope q 
determined by 7' = 7^ through Eq. (|83|) . Note that in this case Qioss = 0. Fig. [3] 
shows the solution at time t = 0.15, obtained with Glimm's method using 
the approximate (left) and exact (right) rarefaction solution in the Riemann 
solver, respectively. Again, from top to bottom, each panel shows gas density, 
velocity, gas pressure and CR pressure. 

Glimm's solutions appear slightly ragged. The raggedness in the right panel is 
purely due to the sampling procedure and does not correspond to an oscillatory 
behavior. In particular, despite their appearance, the points on the rarefaction 
curve remain connected through the correct wave solution and the ragged 
character may or may not be there depending on the specific sample that is 
being drawn. We show an example of this in the next test, where Glimm's 
solution of a shock tube problem is characterized by a smooth rarefaction 
wave. 

On the other hand, the solution on the left panel shows additional irregularity 
which is due to the approximations in Riemann solver. In particular the ap- 
proximations involved in evaluating the solution inside a rarefaction wave add 
spurious structure to the Riemann solution which is picked up by Glimm's 
method. This is undesirable because some of the spurious features may be 
amplified and even become unstable. 

Finally, the left panel of Fig. H] shows the solution obtained with Godunov's 
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Fig. 4. Open symbols: numerical solutions obtained with the hybrid Glimm- 
Godunov's method. Solid line: exact solution. Left panel: rarefaction wave at, 
t = 0.1504683. The initial conditions are described in Eq. (1881 cf. Fig. [3]). Right 
panel: shock tube problem at, t = 0.150794. The initial conditions are described in 
Eq. (I89i). 

method for the same rarefaction wave problem. Overall the numerical solution 
reproduces accurately the evolution of the rarefaction wave. Compared with 
Glimm's method, the solution is now smooth across the wave, although it 
appears slightly less sharp at the head of the wave. 



5.3 Shock Tubes 



We conclude the set of tests with a shock tube problem with the following 
initial conditions 

pi = 1.0, = 0.0, = 1.0, Pi = 0.0, (7^, 

(89) 

pr = 9.0, = 0.0, p; = 10, p; = e.o, = 1.33433. 

The initial CR distribution function for the right state is a power-law, Eq. flHTj) . 
with slope determined by 7^ through Eq. ( !83l) . Note that given the null 
value of the CR pressure on the left state, we need not specify 7^ nor /'(p). 
In the shock tube problem in general a shock, a rarefaction and a contact 
discontinuity develop. In this case the shock is weak and Qioss = 0. 

The numerical results at t = 0.15 are shown in the right panel of Fig. H] 
for the hybrid Glimm-Godunov's method and in Fig. [5] for Glimm's method. 
As in the previous test, the latter figure shows both the results obtained by 
employing an approximate (left) and exact (right) rarefaction solution in the 
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Fig. 5. Opens symbols: shock tube problem solution at, t = 0.150794, obtained 
with Glimm's method, using the approximate (left) and exact (right) rarefaction 
solution in the Riemann solver. Solid lines: 'exact' solution. The initial conditions 
are described in Eq. ([891 cf. Fig. [3|). 

Riemann solver. Note that from top to bottom, each panel now shows gas 
density, velocity, total (gas+CR) pressure and CR pressure, respectively. 

The left panel of Fig. [5] shows additional evidence that when using the approx- 
imate rarefaction solution with Glimm's method, spurious structure appears 
in the numerical solution. This affects not only the rarefaction wave but also 
the intermediate state. When the exact rarefaction solution is employed how- 
ever, Glimm's method produces a highly accurate result. Unlike the earlier 
test, the rarefaction curve is now smooth, confirming how the raggedness in 
the previous test was due to the sampling procedure. Notice how Glimm's 
method preserves the sharpness of both the shock and contact discontinuity. 
As already pointed out, however, their position may be displaced with respect 
to the exact solution by at most one grid cell. 

The right panel of Fig. H] shows that the numerical solution produced by the 
hybrid method is also highly accurate. As a sanity check we notice that in 
switching between Glimm and Godunov's formulations no spurious effect is 
introduced. In addition, the shock position is the same as in Glimm's method 
solution. The rarefaction part of the solution is well captured, although now 
the foot of the rarefaction is not as sharp as in Glimm's method. Finally, 
the contact discontinuity spreads over a few cells, which is characteristic of 
Godunov's method. Note that the contact discontinuity appears both in the 
density and in the pressure components, Pg and Pc- This, however, does not 
affect the total pressure Pg + Pc, which remains perfectly constant in the inter- 
mediate state between the rarefaction and the shock, guaranteeing a correct 
velocity profile as well. 
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6 Conclusions 



In this paper we liave developed a method to include the dynamical effects of 
CR pressure on a hydrodynamical system. 

In smooth flows this is achieved by modification of the characteristic equations 
which define the waves that carry information in the CR-hydro fluid. The 
exchange of energy between the two CR particles and the fluid component as 
a result of diffusive processes both in configuration and momentum space, is 
modeled with a flux conserving method. 

Regarding the solution at shock waves, we have shown that once the accelera- 
tion efficiency has been specified as a function of the upstream conditions and 
shock Mach number, the shock CR mediation and the induced substructure 
can be correctly taken into account in the fluid solution by modifying the 
procedure for the Riemann solver. 

We have implemented two numerical schemes for obtaining time dependent 
solutions on a computational grid. One based on Glimm's method and another 
based on a hybrid Glimm-Godunov method. In both approaches we exploit the 
ability of Glimm's method to preserve the discontinuous character of shocks. 
This is useful because when combined with the aforementioned modified Rie- 
mann solver it provides a natural scheme for advancing the shock solution 
at the correct speed, meaning that the CR dynamical effects are taken into 
account without resolving the shock substructure. 

In smooth flows Glimm's method is not the only possible choice and Go- 
dunov's method can also be employed. Our tests show that Glimm's method 
is rather sensitive to the approximations assumed in the Riemann solver pro- 
cedure. In particular, when evaluating the solution inside a rarefaction wave 
it is important to solve the exact equations (!371) - (!39l) or else spurious features 
will appear. This is not the case with Godunov's method. Compared to the 
version of Glimm's method with the exact rarefaction solution in the Riemann 
solver, Godunov's produces smoother solutions and is only slightly less sharp 
in capturing the head of rarefaction waves. So the hybrid method is also a 
viable option. 

The proposed method can be readily employed for the study of one-dimensional 
models of astrophysical systems, such as simplified radially symmetric de- 
scriptions of, e.g.. Supernova Remnants or Galaxy Clusters. The potential 
benefit of this scheme, however, lies in the possibility of coupling it with 
three-dimensional shock tracking algorithms or extending Glimm's method 
to three dimensions. In this case the dynamical role of CRs in astrophysical 
systems can be studied without geometrical restrictions. Since numerical dis- 
sipation is necessary for stable hydrodynamic shocks in three dimensions 0, 
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but spoils the correct CR-hydrodynamic shock solution, shock tracking algo- 
rithms may offer the only way to self-consistently study CR-hydrodynamics 
in multi-dimensions. 
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